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We present a direct experimental investigation of the thermal ordering in an artificial analogue 
of an asymmetric two dimensional Ising system composed of a rectangular array of nano-fabricated 
magnetostatically interacting islands. During fabrication and below a critical thickness of the mag¬ 
netic material the islands are thermally fluctuating and thus the system is able to explore its phase 
space. Above the critical thickness the islands freeze-in resulting in an arrested thermalized state 
for the array. Determining the magnetic state of the array we demonstrate a genuine artificial two- 
dimensional Ising system which can be analyzed in the context of nearest neighbour interactions. 


The Ising model invented by Wilhelm Lenz and solved in 
one dimension by Ernst Ising in 1924 is one of the pillars 
of statistical mechanics [I, 2]. Although built on a simple 
basis, that of an interacting system composed of a chain 
of entities with only two discrete states, s = {I,—I}, 
the Ising model still to this day is used to model mag¬ 
netic systems and can be applied to a wealth of atomistic 
and mesoscopic experimental systems ranging from ferro¬ 
magnetic ordering and atomic-scale antiferromagnets [3] 
to the ordering of binary colloidal structures [4] and ther¬ 
mal artificial spin systems [5]. In the two dimensional 
case spins are arranged on a square lattice and each spin 
interacts with four nearest neighbours, as seen in Fig. I. 
The total energy of such a two-dimensional Ising system 
can be described by the equation 
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where Jy and Jy correspond to interaction energies in 
the horizontal [10] and vertical [01] lattice directions of 
the two dimensional crystal and the sum is taken over all 
pairs of nearest neighbour spins and Sj. As opposed 
to the one dimensional model which shows no sponta¬ 
neous magnetization at temperatures T > 0 a spon¬ 
taneous magnetization appears in the two dimensional 
case[6] with an order parameter given by[7, 8] 
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where T corresponds to the temperature and is the 
Boltzmann constant. The model is not complicated by 
the choice of different values or signs of Jy and Jy [7]. 
For J > 0 the interaction between neighbouring spins 
favors a parallel alignment (see Fig. 1) while in the case 
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FIG. 1. The two dimensional Ising system. A two dimen¬ 
sional Ising system consists of an array of spins which are 
only capable of pointing in two opposite directions, s = ±1 
and interacting with their four nearest neighbours with in¬ 
teraction energies Jh and Jy along the [10] and [01] lattice 
directions. At T = 0 K the ground state is doubly degenerate 
consisting of all spins pointing in the (+) or the (-) directions 
for the case of a symmetric Ising system. For T > 0 K de¬ 
fects above the ground state occur enclosed by a domain wall 
separating the areas of spins in the all (+) or all (-) directions. 


of J < 0 spins have a preference for an antiferromagnetic 
alignment. 

Nano-patterned single-domain magnetic thin film is¬ 
lands have been a prominent candidate in recent years 
for creating artificial analogues of interacting systems. 
Using modern lithographic techniques it has become fea¬ 
sible to design directly the shape of such islands and 
their geometrical arrangement creating two dimensional 
arrays of interacting artificial spins. Artificial Ising-like 
spins can be realized by designing elongated islands of 
thin films materials, for which shape anisotropy confines 
the magnetization to only two possible orientations. By 
arranging such islands in different geometries a wealth 
of interacting systems can be studied including cellular 
automata [9, 10] and frustrated systems such as artifi¬ 
cial spin ice [11, 12]. The two dimensional nano-scale 
nature of these systems enables their state to be di¬ 
rectly determined by imaging techniques such as mag- 
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FIG. 2. Nano-patterned artificial spins. Atomic force mi¬ 
croscopy (AFM) image of part of the patterned array showing 
the structure and arrangement of the nano-patterned islands 
(see Methods). The shape of the islands confines the magne¬ 
tization to point along the long axis of the islands and dur¬ 
ing thermalization at the onset of a frozen state they can be 
considered as thermally active superspins. In the superspin 
model the vertical coupling ([01] direction) prefers a ferro¬ 
magnetic alignment of the spins (> 0) whereas the lateral 
coupling ([10] direction) prefers an antiferromagnetic align¬ 
ment {Jh < 0) resulting in the ground state ordering of the 
artificial spins illustrated to the right in the figure. 


netic force microscopy [11, 13], photoemission electron 
microscopy [14-16] and Lorentz microscopy [17]. 

In this paper we investigate a two dimensional array 
composed of elongated thin film islands in a square lattice 
(see Fig. 2). During fabrication the array goes through a 
dynamic phase enabling the system to thermally explore 
its phase space leading to a low energy ordered state [13]. 
After this dynamic phase the state of the system becomes 
frozen-in, generating a snapshot of an arrested thermal- 
ized state. The determination of the magnetization of 
individual islands in the array allows us to demonstrate 
that the dynamic phase leads to an ordering of the mag¬ 
netic state of the array, which can be described by the 
two-dimensional Ising model. 

The thickness regime wherein thermal dynamics of the 
spins are obtained occurs during the deposition of the 
magnetic material onto prepatterned substrates as shown 
by Morgan et al [13, 18]. During the deposition, the is¬ 
land thickness (and thereby their volume) becomes grad¬ 
ually larger and the magnetization reversal energy bar¬ 
rier, associated with their shape anisotropy increases. 
In the initial stages of the island growth this barrier is 
smaller than the thermal energy enabling the magnetiza¬ 
tion to spontaneously fluctuate between the two low en¬ 
ergy states defined by the shape anisotropy. As the thick¬ 
ness becomes larger the reversal energy increases, eventu¬ 
ally overcoming the thermal energy, reaching a threshold 
where the superparamagnetic behaviour is suppressed as 
the scale of the combined shape anisotropy reversal en¬ 
ergy barrier and the energy landscape of the array are 
sufficient to lock the magnetization in each of the islands 
(see Fig. 3). Subsequently the array can be imaged by 


magnetic force microscopy (MFM) and the magnetiza¬ 
tion direction of each island can be determined (see Fig. 

4 )- 

During the limited time window where fluctuations oc¬ 
cur the fast dynamics allow the array to explore its phase 
space and achieve an equilibrium condition. During this 
dynamic phase before the magnetization becomes frozen 
the reversal energy barrier and the interaction energies 
between neighbouring islands is of the order of the ther¬ 
mal energy /cpT. For the island sizes and array param¬ 
eters used for this study the energy values involved cor¬ 
respond to room temperature values for film thicknesses 
in the nm regime. Uniform thermal dynamics over 
the entire array therefore require a well defined, stable 
thickness of each island in order to minimize random¬ 
ization effects due to variations in the island thickness 
and film roughness. Amorphous magnetic materials dis¬ 
play soft magnetic properties and a high degree of struc¬ 
tural uniformity and are thereby suitable for well de¬ 
fined layers below 1 nm in thickness [19]. For this study 
we therefore choose amorphous Co 68 Fe 24 Zr 8 as the is¬ 
land material previously used for creating ultra-thin mag¬ 
netic layers [20] and well defined nano-patterned multi¬ 
layered structures [21]. Furthermore, a field imprinted 
anisotropy can be induced in Co 68 Fe 24 Zr 8 enhancing the 
energy barrier for reversal as the magnetization has set¬ 
tled in a fixed direction [22]. 

Considering magnetostatic interactions in the point 
dipole approximation between neighbouring islands re¬ 
veals an interaction scheme which can be mapped to a 
ferromagnetic interaction in the vertical direction Jy > 0 
and antiferromagnetic in the horizontal direction < 0. 
The lowest energy state of the array is therefore com¬ 
posed of a staggered arrangement of ferromagnetically 
aligned chains (see Fig. 2). The lowest energy state of 
an asymmetric two dimensional Ising system is two de¬ 
generate, in our case corresponding to a ferromagnetic 
ordering in the vertical direction and an antiferromag¬ 
netic ordering in the horizontal direction. Excitations 
above the ground state occur through the reversal of a 
macrospin and can be viewed in the form of boundary 
walls separating the two possible ground states as illus¬ 
trated in Fig. 4(c) and the energy state of individual 
island can be categorized into 9 different energy states of 
varying degeneracy shown in Fig. 4(d). 

Counting the abundances of excitations composed only 
of independent vertical or horizontal excitations the rel¬ 
ative energy scale between the two directions can be at¬ 
tained. The observed probabilities of these excitations 
decreases exponentially with increasing energy (inset in 
Fig. 5) in accordance with a Boltzmann distribution of 
the states. Determining the ratio of the excitation en¬ 
ergies from the inset the energies relating to all energy 
values of the individual islands can be listed in units 
of the energy involved with a horizontal excitation \Jh\' 
Within the combined plot (Fig. 5) the observed abun¬ 
dances decrease exponentially establishing the probabil¬ 
ity for a macrospin to be in an energy state E to be 
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FIG. 3. Snapshot of an arrested state. The ordering of the investigated array illustrated through highlighting the two degenerate 
ground state and the domain walls separating them. As compared to a fully ferromagnetic 2-dimensional Ising system the 
antiferromagnetic coupling leads to an additional asymmetry in the population of the up and down islands themselves, which 
is indicative of the influence of an external field. Out of the N — 10487 islands analyzed riu — 5294 were found to point 
in the up direction while rid — 5193 were found to point in the down direction resulting in an asymmetry within statistical 
limits. An overview of the entire investigated array as well as a data hie containing the directions of all islands is supplied in 
Supplementary material. 


given by a Boltzmann distribution ^ exp(F^//cBT) and 
that the system can be sufficiently described by a nearest 
neighbour interaction model. 

Identifying the domain walls separating the two degen¬ 
erate ground states of the array facilitates the mapping 
of the magnetic structure of the system as two ordering 
states, as illustrated in Fig. 3. The mapping of these 
two states onto two domain colours for the entire array 
is shown in Fig. 6. Counting the number of islands falling 
into each of these two domains the order parameter of the 
array can be obtained. The resulting order parameter for 
the array, dehned by M = — riw)/{nb + riw ) where rib 

and riw correspond to the domain populations of black 
and white domains, can then be obtained. Truncating 
the data array to a square shape the array size is reduced 
to n = 9828 out of which rib = 5283 islands fall into the 
black domain while = 4545 islands fall into the white 
domain. The resulting value of M = 0.075 ± 0.015 re¬ 
veals a slightly higher population of the black domains. 
Although the statistics of this value are limited by the 
finite number of observed islands in the array it could 
be indicative of the array being in a state close to, or 
above Tc, since the lack of global order in the macro¬ 
scopic arrangement can lead to a finite value of the order 


parameter. The listed uncertainty of M is determined 
from the square root of the domain populations, and 
rib. 

By calculating the pairwise correlation between spins, 
at a distance r, within the experimental array the cor¬ 
relation function for the system, G(r), can be deter¬ 
mined. The resulting correlation array is shown in Fig. 
7. The preference for an antiferromagnetic arrangement 
of neighbouring spins in the [10] direction introduces the 
possibility of negative values in the correlation and alter¬ 
nating positive and negative values along the [10] direc¬ 
tion. Figure 7 therefore shows the absolute value of the 
pairwise correlation |G(r)| as determined from the array. 
As can be seen in Fig. 7 the pairwise correlation follows 
an exponentially decreasing function as expected for an 
extended array of Ising like spins at T > Tc. Further¬ 
more, the asymmetric nature of the array is revealed in 
the different values of the correlation lengths along the 
[10] and [01] directions of the array with ^[lo] = 2.87 and 
C[oi] = 1.02. The correlation length in the [11] direction, 
C[ii] = 0.97, is revealed to be similar to ^[oi]. 

The role of the ratio of the interaction strengths on 
the order parameter and specifically the ordering tem¬ 
perature of the array can be investigated using Onsager’s 
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FIG. 4. Macrospin arrangements imaged by MFM. Analysis 
of a part of the MFM images recorded showing a portion of 
the N = 10487 island array for which the direction of the mag¬ 
netization was determined. The lateral extent of the entire 
nano-patterned array was 2x2 mm^, corresponding to a to¬ 
tal number of islands 40 million, a Magnetic force microscopy 
image showing part of the analyzed array, b Schematic high¬ 
lighting of island contours and the excitation boundaries, c 
The energy state of the islands quantified into the nine differ¬ 
ent possible energy states of individual islands with respect 
to the orientation of their nearest neighbours shown in d. 
The energy states along with their respective degeneracies are 
listed with respect to the lowest energy ground state assuming 
that Jh JV' 


solution, as initially done by Chang [23] or utilizing nu¬ 
merical calculations such as Monte Carlo simulations, see 
Fig. 8. In the case of an isotropic system, wherein the in¬ 
teraction strength between neighbours, J, is the same in 
the two main lattice directions ([01] and [10]), the order¬ 
ing temperature is given by Tc = « 2.269J//cb. 

As Jy decreases with respect to Jh the relative Ty of 



FIG. 5. Macrospin energy state statistics. The number of 
observations of each of the 9 energy states illustrated in Fig. 
4(d) divided by their degeneracy. The inset shows the num¬ 
ber of observed energy states composed exclusively of either 
vertical excitations (energy levels 2Jv and 4J^;), blue line, or 
horizontal excitations (energy levels 2Jh and 47/^), red line. 
From the ratio of the slopes the ratio of the energy scale 
between the vertical and horizontal interaction energies is de¬ 
termined to be \Jv/Jh\ = 0.2943. Gonsidering this ratio the 
number of observations of the 9 different energy states can be 
listed as a function of their individual energy in units of 2| 

The error bars correspond to the square root of the number 
of observations for each of the energy states. 



FIG. 6. Macrospin array domain configuration. The 
macrospin conhguration of the investigated arrays mapped to 
black and white domains corresponding to the two degenerate 
ground state of the system. 


the system also decreases. Figure 8 shows results from 
Monte Carlo simulations for decreasing values of the ra¬ 
tio \Jv/Jh\‘ The change in Ty corresponds well to what 
is expected in the two dimensional Ising model as the 
probability of observing two parallel spins in the [01] di¬ 
rection decreases. As can be seen in Fig. 8 the results 
of numerical simulations correspond well to the analyti¬ 
cal solution with small deviations arising from the finite 
nature of the simulated array. 

Assuming a system in a perfectly ordered ground state 
at 0 K as the temperature is increased thermal energy is 
introduced into the system and spins will start to change 
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FIG. 7. Pairwise correlations G for the array. The graph 
shows the absolute value, |G|, for easier mapping of the pair¬ 
wise correlations for both the ferromagnetic and antiferro¬ 
magnetic directions. From the slopes the resulting correla¬ 
tion lengths along the different directions is determined to be 
^[ 10 ] = 2.87, ^[ 01 ] = 1.02, and ^[n] = 0.97 




T=1.5Jh/kB 



"'TT2.0Jh/kB' 


their directions. The perfect initial order is therefore bro¬ 
ken and the system moves towards a thermally disordered 
state with a reduced order parameter. From numerical 
simulations this transition can, furthermore, be observed 
through the pairwise correlation G(r) between spins. For 
the two dimensional Ising model the correlation function 
can be written as G(r) ^ exp( —|) for T < Tc and T > Tc 
where ^ is the correlation length affected by temperature 
and the interaction strength between neighbouring spins. 
The inset in Fig. 8 shows the temperature dependence 
of ^ for the three major directions of the array, [10], [01], 
and [11] for the asymmetry ratio \Jv/Jh\ = 0-3 as deter¬ 
mined from Monte Carlo calculations. As Jh > Jy the 
corresponding correlation length in the [10] direction is 
larger than in the other directions while the correlation 
length in the [01] and [11] directions are similar. Con¬ 
sidering the results of the Monte Carlo simulations for 
an interaction ratio of \Jv/Jh\ = 0.3 one can see that 
for temperatures larger than Tc the correlation lengths 
^[ 10 ] and ^[ 01 ] follow the same trend as the experimental 
results, i.e. C[io] > ^[oi] and ^[oij r\j ^[ 11 ]. In particular 
at T = 1.9Jh/kB the values obtained from the simula¬ 
tions are ^[lo] = 1.95, and [01], ^[oi] = 1.06, ^[n] = 1.08, 
similar to the values obtained experimentally. These re¬ 
sults indicate that in the arrested state the experimental 
observations can be described with an effective nearest 
neighbour Ising model. Hence, long range dipolar inter¬ 
actions are not needed to model the present system in the 
arrested stated. However, dipolar interactions cannot be 


FIG. 8. Order parameter versus asymmetry, a Results of 
Monte Garlo simulations showing the order parameter M as 
a function of temperature, in units of Jh/^B, for different val¬ 
ues of \Jv/Jh\- Each point in the graph corresponds to an 
average of several realizations (see methods). As the ratio 
decreases from the symmetric value of the \Jv/Jh\ = 1 the 
Tc also decreases. The analytical Onsager solution given by 
equation (2) (solid lines) fits quite well with the obtained 
numerical data. The inset shows the correlation length ^ 
as a function of temperature for the ratio \Jv/Jh\ — 0.3 
for the three major crystallographic directions of the array, 
[10], [01], and [11]. b - e Representative snapshots of the nu¬ 
merical simulations at temperatures of 1.0 J/i//cb, 1.34J/i//cb, 
1.5J/i//cb, and 2.0 J/i//cb, respectively, for an asymmetry ra¬ 
tio of \Jv/Jh\ — 0.3. Due to the asymmetry the correlation 
length in the [10] direction is larger than in the [01] direction 
stretching out the domains in the [10] direction, similar to 
what is observed in the experimental data. 


completely excluded if one desires a full description of 
the system in any non-arrested state. 

Further advances in materials science and magnetic imag¬ 
ining techniques will undoubtedly allow us to follow 
the thermal evolution of a multitude of different ar¬ 
tificial spin systems[24] such as artificial spin ice, one 
and two dimensional Ising arrays and different frustrated 
arrangements[25, 26] as they explore their phase space. 
Such investigations do not merely offer a direct real space 
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probe of the microstate of thermal systems at a local scale 
but also a direct determination of the dynamics involved. 
The possibility of determining directly the state of these 
systems enables the experimental probing of the effect of 
external variables such as temperature and applied mag¬ 
netic field. This introduces the possibility of investigating 
e.g. the two-dimensional Ising model under applied ex¬ 
ternal field directly for which an analytical solution does 
not exist. 


Methods 

The island array was prepared by magnetron sputtering 
thin film growth onto a pre-patterned fused silica sub¬ 
strate. The 2x2 mm^ pre-patterned array was fabricated 
by surface conformal nano-imprint lithography [27, 28]. 
The array was composed of 470 nm x 160 nm islands 
arranged in a square lattice with a periodicity of 200 nm 
along the short axis of the islands ([10] direction) and a 
periodicity of 500 nm along the [10] direction. The amor¬ 
phous film was deposited by magnetron sputtering from a 
compound target at room temperature. Argon was used 
as the sputtering gas at a pressure of 3.0 mTorr. Be¬ 
fore growth the base pressure of the sputtering system 
was below 2 x 10“^ Pa. Initially a 2 nm thick layer of 
AlyoZrso was grown onto the pre-patterned substrate fol¬ 
lowed by the growth of a 7 nm thick Co 68 Fe 24 Zr 8 and 
finally a 2 nm thick Al 7 oZr 3 o layer to prevent degrada¬ 
tion. After patterning the structural quality of the films 
was investigated by AFM carried out in contact mode 
using a Nanosurf Mobile S instrument. 

The thermally ordered state of the array was investigated 
by magnetic force microscopy using a digital instrument 
dimension 3100. Topography was imaged by height con¬ 
trast in tapping mode and magnetic micrographs were 
scanned in lift mode with a lift scan height of 70 nm. 
In order to reduce the risk of the stray field from the 
MFM tip altering individual island states the images were 


recorded using a low magnetic moment Co alloy MFM tip 
(PPP-LM-MFMR). Repeated scanning of the same area 
did not alter the state (see supplementary material). 
The system was modeled as an Ising system, in which ev¬ 
ery nano-island is considered as a single macrospin with 
only two possible orientations (+1 and —1), the simula¬ 
tion cell consists of 150 x 150 spins. Simulations with both 
periodic and open boundary conditions were performed. 
The interactions were only considered between nearest 
neighbours. The properties of the system were stud¬ 
ied via Monte Carlo simulations using the Metropolis- 
Bastings algorithm. To simulate the process of the spins 
freezing at a given configuration and then evolving to¬ 
wards the most stable configuration before any measure¬ 
ment is performed the system is considered to be com¬ 
pletely disordered at a temperature higher than its crit¬ 
ical temperature. Tc. During this thermalization phase 
no measurements are performed, the system is allowed 
to evolve in accordance to the usual Metropolis-Bastings 
algorithm. After a certain number of Monte Carlo steps 
the system is cooled down and the process is repeated 
once more. This procedure is then repeated until the 
desired measurement temperature is reached. 
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